#用McCabe方法进行膜级联设计
#根据给定的VRR值进行设计
import numpy as np 
import matplotlib.pyplot as plt

VRR = 3.19#F/R
Rej = np.array([0.99,0.12])
RetA = Rej[0]
RetB = Rej[1]
Abta = 1-1/(1-1/VRR)*(1-(1/VRR)**RetA)
Abtb = 1-1/(1-1/VRR)*(1-(1/VRR)**RetB)
alpha = (1-Abtb)/(1-Abta)*(1-(1-1/VRR)*(1-Abta))/(1-(1-1/VRR)*(1-Abtb))

Cin = np.array([10,2])#进料的A/B（染料/盐）组分浓度
xf = Cin[1]/np.sum(Cin)#盐比例
xtarget = 0.96
xin = 1-xtarget

xps = []#渗透侧集合
xrs = []#渗余侧集合

count = 1#总膜数计数器
countper = 0#渗透段计数器
countret = 0#渗余段计数器

xi = xin
yi = xi 
Cper = Cin
C = Cper#迭代初值
Gamas = []
Ins =[]

plt.figure(1)
plt.plot([0,1],[0,1],'k')#对角线
plt.xlabel('Retentate mole fraction(xR)')
plt.ylabel('Permeate mole fraction(xp)')
#膜分离线
def equil(x):
	y = alpha*x/(1-x)/(1+alpha*x/(1-x))
	return y
while True:
	if xtarget <= xi:
		break
	xt = xi
	Theta = (1-1/VRR)*xi*((1-Abtb)+(1/xi-1)*(1-Abta))
	Gama = (1-Theta)/Theta 
	Gamas = np.append(Gamas,Gama)
	T = yi/(1-yi)/alpha 
	xi = T/(T+1)#平衡线的x
	plt.plot([xt,xi],[yi,yi],'b')#横线
	xrs = np.append(xrs,xi)
	In = 0
	for i in range(1,count+1):
		G = Gamas[0:i]
		In += np.sum(np.cumprod(G))
	Ins = np.append(Ins,In)
	#渗透侧操作线方程组
	def per(x,I):
		y = I/(I+1)*x+1/(1+I)*xin 
		return y 
	#渗余侧操作线方程组
	def ret(x,I):
		y = (1+1/I)*x-1/I*xtarget
		return y
	xps = np.append(xps,yi)
	I = In
	xinter = I*(1+I)/(2*I+1)*(xin/(1+I)+xtarget/I)#本级操作线交点
	plt.plot([xin,xinter],[xin,per(xinter,I)],'r')#渗余侧操作线
	plt.plot([xtarget,xinter],[xtarget,per(xinter,I)],'g')#渗透测操作线
	#判定属于哪一段操作
	if xi <= xinter:
		yi = per(xi,I)
		countper += 1
	else:
		I = Ins[count-1]
		countret += 1
		yi = ret(xi,I)
	plt.plot([xi,xi],[equil(xi),yi],'b')#竖线
	count += 1
xeq = np.arange(0,1,0.001)
yeq = []

for j in xeq:
	yeq = np.append(yeq,equil(j))

plt.plot(xeq,yeq,'y')
plt.plot([xin,xtarget,xf],[xin,xtarget,xf],'o')#三个点
plt.xlim((0, 1))
plt.ylim((0, 1))
print([count,countper,countret])
print(xps)
print(xrs)
plt.show()


	

	
	